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Abstract-  The  Davies-Bouldin  cluster  validation  technique  has 
been  utilized  to  determine  the  optimal  time  series  parameters 
(model,  order  and  number  of  coefficients)  to  be  used  in  a 
myoelectric  pattern  recognition  system  based  on  intramuscular 
recordings. 

The  data  used  in  our  longitudinal  study  were  recorded 
intramuscularly  (4  muscle  channels)  from  an  amputee  subject 
over  three  different  experimental  days.  During  the  recording 
sessions,  the  subject  performed  4  wrist  movements  (pronation, 
supination,  wrist  flexion,  wrist  extension)  in  a  way  that  was 
intuitive  to  her.  The  time  series  coefficients  obtained  from  the 
EMG  data  (from  AR,  MA,  and  ARMA  modeling)  were  used  for 
generating  feature  spaces  with  4  classes  (1  per  movement).  The 
results  showed  that  the  optimal  parameters  didn’t  differ 
substantially  (AR,  order  4,  2  coefficients)  for  those  obtained  for 
other  type  of  movements  and  surface  EMG  recordings.  The  use 
of  the  Davies-Bouldin  technique  during  the  analysis  produced 
important  information  about  the  separability  and  consistency  of 
the  data  across  different  days. 

Keywords  -  Myoelectric  control,  EMG,  intramuscular 
electrodes,  pattern  recognition,  time  series,  cluster  validation, 
Davies-Bouldin  index. 

I.  Introduction 

Multifunctional  myoelectric  control  usually  entails 
cumbersome  operation  and  time  consuming  training. 
Throughout  the  years,  different  pattern  recognition 
techniques  and  EMG  feature  spaces  have  been  evaluated 
with  the  purpose  of  solving  the  two  above-mentioned 
problems.  Some  of  the  most  widely  used  features  have  been 
based  on  time  series  coefficients  derived  from  autoregressive 
(AR),  moving  average  (MA),  and/or  autoregressive-moving 
average  (ARMA)  models  [2,3].  The  rationale  for  the  use  of 
time  series  coefficients  is  based  on  the  relationship  between 
coefficients  and  the  EMG  frequency  content.  Likewise, 
frequency  content  is  movement  dependent  due  to  changes  on 
the  motor  units  populations  involved  in  performing  different 
movements.  Because  of  computational  convenience  [2],  the 
preferred  model  type  and  number  of  coefficients  utilized 
with  surface  EMG  has  been  AR  with  order  3  or  4. 
Nonetheless,  some  investigators  [1]  have  claimed  that  time- 
series  modeling  of  EMG  signals  should  not  be  considered 
static  or  invariant  because  the  spectral  behavior  of  EMG  data 
is  dependent  on  the  specific  muscle,  contraction  level,  and 
limb  function.  However,  their  efforts  revealed  that  the  best 
model  representing  surface  EMG  was  still  AR  with  order  4. 
Because  of  the  differences  in  our  study  (wrist  movements, 
intramuscular  EMG,  EMG  feature  extraction  performed  at 
the  contraction’s  onset  as  opposed  to  extraction  during 
sustained  contraction),  we  believe  that  the  results  from 
former  studies  may  not  be  applicable  to  our  problem.  Thus, 


we  performed  an  analysis  to  identify  the  optimal  time-series 
parameters.  With  this  purpose  we  utilized  the  Davies- 
Bouldin  cluster  validation  technique  [4].  This  crisp  cluster 
validation  technique  presents  highly  desirable  characteristics 
to  solve  the  problem  at  hand:  (1)  it  performs  an  evaluation 
based  on  the  geometrical  characteristics  of  the  feature  space 
(it  identifies  compact  and  well-separated  partitions  of  the 
feature  space),  and  (2)  it  has  reduced  sensitivity  to  cluster 
outliers  [5].  In  addition,  other  authors  [6,  7]  have 
successfully  applied  the  Davies-Bouldin  index  to  perform 
comparisons  between  feature  spaces  based  on  EMG  data. 

II.  Methodology 

A.  Experimental  protocol  and  pre-processing 

The  subject  was  a  female  (age  33)  presenting  a  below- 
elbow  amputation  of  her  left  arm  (non-dominant  hand)  and 
strong  phantom  limb  sensation.  The  subject  never  used  a 
prosthesis  since  the  time  of  the  amputation  (13  years  before 
the  experiments). 

A  pair  of  coiled  wire  electrodes  for  bipolar  recording  was 
inserted  with  2cm  tip  separation  into  each  of  the  following 
muscles:  (1)  Pronator  Teres,  (2)  Supinator,  (3)  Flexor  Carpi 
Radialis,  (4)  Extensor  Carpi  Radialis.  During  a  period  of  30 
days,  the  electrodes  remained  chronically  implanted  and, 
within  that  period,  recording  sessions  were  performed  on  12 
different  days.  In  our  study  we  used  data  obtained  during  the 
last  3  experimental  sessions  (i.e.,  the  10th,  11th,  12th  sessions). 
The  subject  was  seated  facing  a  computer  display  that 
executed  the  animations  of  the  movements  to  be  evaluated 
(Fig.l).  For  each  given  movement,  the  subject  was  requested 
to  follow  the  animations  of  the  hand  while  contracting  her 
residual  muscles  during  4  blocks  of  exercises  (25  trials  each). 
The  EMG  signals  were  amplified  (between  1000  and  10000 
times)  and  bandpass  filtered  (10Hz  to  1kHz)  before  being 
sampled  at  2kHz.  This  sample  rate  was  selected  due  to  the 
bandwidth  of  intramuscular  recordings  (1kHz). 

For  computational  efficiency,  we  segmented  the  data  into 
separated  files  that  contained  at  least  4  trials.  We  ran  an 
onset-detection  algorithm  on  the  trials  contained  in  these 
files  until  we  obtained  the  following:  (session  10)  51  onsets 
per  movement,  (session  11)  52  onsets  per  movement, 
(session  12)  50  onsets  per  movement,  which  we  considered 
sufficient  to  perform  our  analysis.  The  algorithm  applied  to 
perform  onset  detection  was  based  on  a  thresholding 
technique  that  evaluated  the  spatial  as  well  as  temporal 
characteristics  of  the  EMG  data  to  perform  a  decision  (see 
[8]  for  more  details).  For  each  given  movement,  after  an 
onset  was  detected,  a  segment  of  200ms  was  collected  from 
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each  muscle  channel  and  successively  stored  in  an  “onset” 
matrix  (4  columns,  1  per  muscle  channel). 
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Fig.  1  Experimental  setup 


B.  Time  series  analysis  using  Davies-Bouldin  cluster 
validation  technique 


(1)  Individual  analysis:  Initially,  we  performed  an 

individual  analysis  on  the  data  obtained  from  each 
experimental  day.  For  a  given  wrist  movement,  and 
following  each  onset  detection,  we  calculated  the  time  series 
coefficients  for  each  of  the  muscles.  These  were  the 
columns  in  an  “onset”  matrix.  We  generated  the  feature 
spaces  separately  for  the  AR,  MA,  and  ARMA  models.  For 
each  model  type  we  repeated  the  feature  extraction  while 
uniformly  increasing  the  order  and  number  of  coefficients. 
The  Davies  Bouldin  indices  (DBIs)  were  calculated  for  each 
of  the  resulting  feature  spaces,  as  follows: 
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where  Rn  is  the  cluster  to  cluster  similarity,  S  the  dispersion 
of  the  ith  cluster.  My  is  the  Minkowski  distance  between  the 
centroids  of  clusters  ith  and  jth  and  Ri=max(Rij)  i^j.  In  our 
analysis  we  utilized  q=2  and  p=2.  Thus,  S,  became  the 


Euclidean  distance  between  the  points  of  a  cluster  to  the 
centroid  of  the  cluster  and  My  the  Euclidean  distance 
between  centroids. 

For  each  given  model  and  each  experimental  day,  the 
resulting  DBIs  were  stored  on  matrices  sized  10x10  (rows= 
increasing  order  of  model;  columns=increasing  number  of 
coefficients).  The  limit  of  order  10  for  each  model  type  was 
chosen  to  maintain  the  evaluation  within  the  working  limits 
of  modern  DSP  processors.  In  addition,  the  use  of  up  to  10 
coefficients  allowed  us  to  evaluate  trends  in  the  results.  In 
the  case  of  the  ARMA  model,  the  total  number  of 
coefficients  was  20  in  order  to  allow  a  uniform  contribution 
of  both  types  of  coefficients  and  to  obtain  results  that  could 
be  compared  with  those  obtained  with  the  AR  and  MA 
models. 

(2)  Combined  analysis:  In  order  to  determine  which 
configuration  behaved  in  an  optimal  way  during  the  three 
experimental  days,  we  pooled  all  the  onsets  utilized  during 
the  individual  analysis  in  a  single  set.  Subsequently,  we 
proceeded  to  evaluate  the  feature  spaces  for  each  given 
model  and  a  uniformly  increasing  model  order  and  number 
of  coefficients.  Again,  the  resulting  DBIs  were  stored  in 
10x10  matrices. 

We  were  able  to  obtain  only  minimal  improvements  on 
class  separability  following  significant  increases  on  the  order 
of  model  and  number  of  coefficients.  Because  larger  orders 
and  number  of  coefficients  are  computationally  more 
expensive,  we  used  a  tolerance  value  to  determine  which  was 
the  optimal  configuration  from  a  computational  as  well  as  a 
class  separation  viewpoint.  This  tolerance  value  was 
calculated  from  each  DBI  matrix  as  the  increase  of  5%  over 
the  minimum  DBI  index  (i.e.,  tolerance  value  =  1.05  x 
(minimum  DBI)).  We  chose  the  value  of  5%  because  we  felt 
that  an  increase  in  5%  over  the  minimum  DBI  will  not 
greatly  degrade  class  separability,  but  it  will  significantly 
improve  the  performance  of  the  resulting  system.  Thus,  for 
each  given  DBI  matrix,  once  the  minimum  DBI  value  was 
found,  we  used  the  tolerance  to  find  the  configuration  with 
the  smallest  order  and  number  of  coefficients  that  had  a  DBI 
below  the  tolerance  value. 

III.  Results 

A.  Individual  analysis 

The  smallest  DBIs  were  continuously  obtained  with  AR 
models  (Table  I).  In  this  case,  the  indices  ranged  from  1.155 
to  3.304.  In  the  case  of  MA  and  ARMA,  the  indexes  ranged 
from  6.631  to  23.678,  and  from  7.206  to  23.951, 
respectively.  From  these  values  we  can  observe  that  the 
DBIs  obtained  by  the  feature  spaces  based  on  AR  model 
didn’t  overlap  with  those  obtained  by  MA  and  ARMA 
models. 

A  more  detailed  study  of  the  indices  revealed  that  above 
a  certain  order  and  number  of  coefficients,  the  feature  spaces 
with  a  larger  number  of  coefficients  also  presented  larger 
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DBI  values  (Fig.  2).  This  behavior  was  observed  for  orders 
larger  than  4  and  number  of  coefficients  larger  than  1.  In 
addition,  we  noted  that  for  increasing  orders  of  the  model, 
the  DBI  indices  attained  larger  values  (Fig.  3).  This  behavior 
presented  small  fluctuations,  but  it  become  more  and  more 
evident  as  a  larger  number  of  coefficients  was  utilized.  This 
behavior  was  observed,  depending  on  the  experimental  day, 
for  orders  larger  than  3,  4,  or  5,  and  number  of  coefficients 
larger  than  or  equal  to  1  and  larger  than  2. 

Finally,  when  we  applied  the  tolerance  criteria  to  the  AR 
model,  we  obtained  the  optimal  configuration  as  follows: 

day  10:  order  4  and  2  coefficients  (DBI=1.210). 
day  11:  order  4  and  2  coefficients  (DBI=1.785). 
day  12:  order  4  and  2  coefficients  (DBI=1.512). 


Table i 

Wrist  movements:  means  and  standard  deviations  of  the  DBI  values 
obtained  in  each  experimental  day  for  each  of  the  models:  AR,  MA;  ARMA. 


daylO 

dayl  1 

dayl  2 

AR 

mean 

1.322 

1.951 

1.68 

std 

0.21 

0.225 

0.277 

MA 

mean 

10.544 

11.678 

12.21 

std 

2.65 

2.735 

3.063 

ARMA 

mean 

10.19 

11.343 

10.507 

std 

2.348 

2.606 

1.995 

AR  model  session  10:  wrist  movements 
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Fig. 2  Experimental  day  10:  Dependence  of  the  DBIs  with  the  number  of 
coefficients. 


Fig. 3  Experimental  day  10:  Dependence  of  the  DBIs  with  the  order  of 
model.  Note,  the  configurations  with  DBI  slightly  smaller  that  the  one 
presented  by  the  optimal  configuration  (order  2,  4  coefficients). 

B.  Combined  analysis 

Once  again  the  feature  spaces  based  on  the  AR  model 
produced  the  smallest  DBI  values  (Table  II).  This  time  the 
values  ranged  from  1.6297  to  3.399.  Note,  that  the  spread  of 
the  data  and  the  mean  and  standard  deviations  don’t  differ  to 
a  large  extend  from  the  results  obtained  during  the  individual 
analysis.  The  range  of  values  derived  from  the  MA  and 
ARMA  models  was  from  13.844  to  31.797  and  14.283  to 
31.139,  respectively. 


Table  ii 

Wrist  movements:  means  and  standard  deviations  of  the  DBI  values 
obtained  from  the  pooled  coefficients  obtained  during  all  the  experimental 
days.  Results  obtained  for  AR,  MA,  and  ARMA  models. 


AR 

MA 

ARMA 

wrist 

mean 

1.851 

18.081 

17.979 

std 

0.343 

4.325 

3.901 

We  noted  once  more  that  the  DBIs  obtained  from  the 
feature  spaces  based  on  the  AR  model  didn’t  overlap  with 
those  derived  from  the  MA  or  ARMA  models.  Subsequently, 
we  applied  the  tolerance  criteria  to  the  AR  model  and  found 
that  the  optimal  DBM  value  was  order  4  and  2  coefficients 
(DBM=1.662).  Note  that  this  configuration  is  the  same  for 
all  the  experimental  days  during  the  individual  analyses. 

Furthermore,  we  observed  the  same  relationship  between 
the  number  of  coefficients  and  the  DBIs  that  was  already 
observed  during  the  individual  analysis.  In  this  case,  the 
relationship  appeared  for  orders  larger  than  4  and  number  of 
coefficients  larger  than  1.  We  also  noticed  that  when  the 
order  of  the  model  increased,  the  DBIs  increased  as  well. 
This  behavior  was  observed  for  orders  larger  than  3  and  a 
number  of  coefficients  larger  than  2. 

IV.  Discussion 

The  use  of  an  AR  model  has  been  favored  by  the  results 
from  both  individual  and  combined  analyses.  Additionally, 
the  results  (model  type,  order,  and  number  of  coefficients) 
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were  consistent  with  those  obtained  by  other  investigators 
using  surface  EMG  data  and  more  standard  methods  of  time- 
series  identification  [1,  2],  The  reasons  for  this  likeness 
between  the  results  obtained  with  surface  EMG  recordings 
and  intramuscular  EMG  might  be  due  to  the  low  recruitment 
level  produced  at  the  onset  of  the  muscle  contraction. 
Because  the  units  with  lower  conduction  velocity  (lower 
frequency  components)  are  recruited  first,  the  spectra  of 
intramuscular  recordings  at  the  onset  of  the  contraction 
might  be  similar  to  the  spectra  of  surface  EMG  recordings 
during  sustained  contraction. 

We  also  noted  that  the  model  type,  order,  and  number  of 
coefficients  that  provided  an  optimal  class  separation  was  the 
same  for  the  individual  and  combined  analysis:  order  4  and  2 
coefficients.  Additionally,  in  both  cases  (individual  and 
combined  analysis)  the  spread  of  the  DBI  indices  and  the 
mean  and  standard  deviation  values  were  comparable.  These 
facts  indicate  that  the  distribution  of  samples  in  the  different 
feature  spaces  was  rather  consistent  across  different  days. 
This  result  is  important  when  we  consider  the  use  of  this 
feature  space  in  a  pattern  recognition  system  since,  in  this 
case,  it  would  be  possible  to  obtain  a  consistent  performance 
of  the  system  across  different  days. 

V.  Conclusion 

The  Davies-Bouldin  cluster  validation  technique  has 
permitted  us  to  determine  the  optimal  configuration  (model 
type,  order,  and  number  of  coefficients)  from  a  pattern 
recognition  standpoint.  We  have  obtained  results 
comparable  to  those  from  surface  EMG  studies  although  we 
recorded  intramuscularly,  and  the  type  of  movements  and 
muscles  were  different  from  previous  studies.  The  analysis 
of  DBIs  from  the  different  feature  spaces  generated 
important  information  about  the  separability  and  consistency 
of  the  data  across  different  experimental  days.  Also,  the  DBI 
values  permitted  us  to  evaluate  the  comparative  separability 
of  AR  versus  the  MA,  and  ARMA  models.  Best  results 
overall  were  obtained  with  an  AR  model  of  order  4  and  2 
coefficients. 
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